AC-AliJ  792  AN  0(N(3)L)  PRIHAL-DUAL  INTERIM  POINT  AL60RITHH  PM 
LINEAR  PROGRA«hING<U>  CALIFORNIA  UNIV  BERKELEV 
OPERATIONS  RESEARCH  CENTER  R  C  HONTEIRO  ET  AL.  NAV  07 
UNCLASSIFIED  0RC-S7-4  NBM14-87-K-B2B2  F/0  1274 


ymm 


AD- A 183  792 


DHC.£|l£  COW 


AN  0(n\)  PRIMAL-DUAL  INTERIOR  POINT 
ALGORITHM  FOR  LINEAR  PROGRAMMING** 

by 

MM 

R.  C.  Monteiro 
I.  Adler**** 


ORC  87-4 


May  1987 


This  research  was  sponsored  in  part  by  the  United  States  Navy  Office  of 
Naval  Research  under  contract  N00014-87-K-0202. 


Department  of  Industrial  Engineering  and  Operations  Research.  University  of 
California  at  Berkeley.  Berkeley,  California,  94720. 


L*  ..  .  ,-t  f.  4  i 

-  ■  **,  jj 


87 


8 


i  10  7 


MCWMTV  CLASSiriCATlOM  Of  THIS 


REPORT  DOCUMEMTATtON  RACE 


had  nmucnom 

BEFORE  COMFLtTOO  FORM 


OPC  87-4 


«.  Tint  fan*  AvAmfaJ 

3 

An  0(n  L)  Primal-Dual  Interior  Point  Algorithm 


for  Linear  Programming 


■  AyTHORfaJ 

R.  C.  Honteiro 
I.  Adler 


■  FERCORMINO  OROANIIATION  NAMt  AND  A 

Operations  Research  Center 

3115  Etcheverry  Hall  - 

U.  C.  Berkeley  Berkeley,  CA  94720 


n.  con  troll  in  0  or  net  namk  and  adorcu 

Office  of  Naval  Research 

Dept,  of  the  Navy 
Arlington,  VA  2P2J7 


Technical 


FBRrORHINO  OR*.  RtRORT  NUMBER 

87-4 


I.TniTH  14 


N000  14-87-K-0202 


It.  RCRORT  DATA 

May  1987 


CmimIIAM  OHlM> 


curitt  clam.  r*r  mu  n*«v 


igummmma 


It  DISTRIBUTION  STATEMENT  (ml  mu  R  apart; 


Distribution  Unlimited;  Approved  for  public  release. 


IT.  DISTRIBUTION  STATEMENT  (ml  m»  mturmti  mini<N  WmI  If,  w  EIInmI  *«b  R«mD 


IS.  NET  WORDS  (Cmntlnum  an  rmrmitm  mltm  U  nmmmtmmry  «f  WmHIT  Af  tt—k 


SO  ABSTRACT  (CwiIIdm  an  rnraraa  NR  IlMiuHq  *f  Itmmutf  tf  UmI 

(see  report) 


I  JAN  71  1473  coition  or  I  NOV  ••  11  oobolbti 

»/N  0102-  LA- 014- 1101 


MCwRlTV  CLAMtClC 


An  0(n3L)  Primal -Dual  Interior  Point 
Algorithm  for  Linear  Programming 


R.  C.  Monteiro  and  I.  Adler 


Department  of  Industrial  Engineering  and  Operations  Research 
University  of  California 
Berkeley,  CA  94720 


n 


t  <sv, 


May  1987 

Abstract  -  describe  a  primal-dual  interior 
gramming  problems  which  requires  a  total  of  O(n^L)  arithmetic  operations.  Each 
iteration  updates  a  penalty  parameter  and  finds  an  approximate  .Newton’s  direc- 
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tion  associated  with  the  Kuhn-Tucker  system  of  equations  whi^fi  characterizes  a 
solution  of  the  logarithm  barrier  function  problem.  This  direction  is  then  used  to 
find  the  next  iterate.  The  algorithm  is  based  on  the  path  following  idea.  The 


total  number  of  iterations  is  shown  to  be  of  the  order  of 1 


A 
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Key  Words  -^  Interior-point  methods,  Linear  Programming,  Karmarkar’s  LP 
algorithm.  Polynomial-time  algorithms.  Barrier  function,  Path  following. 


1.  Introduction 


Consider  the  linear  programming  problem 


(P)  min  cTx 
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where  A  €  Rmx".  This  paper  presents  an  algorithm  for  Linear  Programming  (LP)  problems 
based  on  the  logarithmic  barrier  function  approach.  The  logarithmic  barrier  function 
method  was  first  used  for  LP  problems  by  Frisch  [1J.  The  introduction  of  the  new  interior 
point  algorithm  by  Karmarkar  [4]  led  researchers  to  reconsider  the  application  of  the  loga¬ 
rithmic  barrier  function  method  to  LP  problems.  Recently,  this  method  was  first  considered 
by  Gill  et  al.  (2]  to  develop  a  projected  Newton  barrier  method  for  solving  LP  problems. 
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Meggido  [6]  provides  a  theoretical  analysis  for  the  logarithmic  barrier  method  and  proposes 
a  primal-dual  framework  based  on  a  consideration  of  a  pair  of  primal  and  dual  LP  prob¬ 
lems.  Kojima  et  al.  [S],  using  this  framework,  present  an  algorithm  that  works  simultane¬ 
ously  with  a  pair  of  primal  and  dual  LP  problems.  Their  algorithm  is  shown  to  converge  in 
at  most  O(nL)  iterations  with  a  computational  effort  of  0(ft3)  arithmetic  operations  per 
iteration,  resulting  in  a  total  of  0(rt*L)  arithmetic  operations. 

In  this  paper,  we  build  on  the  ideas  in  [S]  and  obtain  a  faster  algorithm.  The  primal- 
dual  framework  presented  in  [6]  is  used.  The  directions  generated  by  our  algorithm  are 
essentially  the  same  as  the  directions  generated  by  the  algorithm  of  Kojima  et  al.  [5].  How¬ 
ever,  working  closer  to  the  "path  of  solutions"  (c.f.  [5]),  we  are  able  to  obtain  convergence 

in  at  most  O(^fnL)  iterations.  Each  iteration  involves  the  inversion  of  a  n  x  n  matrix 
which  can  be  done  in  at  most  0(n3)  arithmetic  operations.  Based  on  ideas  presented  by 
Gonzaga  [3],  we  are  able  to  exploit  the  special  structure  of  the  matrix  to  be  inverted  so  that 
it  can  be  done  in  an  average  of  0(n2-5)  arithmetic  operations  per  iteration.  Thus  overall  our 
algorithm  requires  0(n3L)  arithmetic  operations.  It  should  be  noted  that  Renegar  [7]  was 

the  first  to  introduce  an  interior  point  algorithm  requiring  0(*JnL)  iterations  and  0(n35L) 
arithmetic  operations.  Subsequently,  Vaidya  [8]  improved  it  so  that  the  total  complexity  is 
0(n3L)  arithmetic  operations.  Equivalent  complexity  was  also  obtained  by  an  algorithm 
which  was  presented  by  Gonzaga  [3].  Both  Vaidya’s  and  Gonzaga’s  algorithms  are  primal 
algorithms.  It  should  be  noted  that  in  order  to  simplify  the  complexity  analysis  presenta¬ 
tion,  we  assume  throughout  the  paper  that  m  =  O(n). 

Our  paper  is  organized  as  follows.  In  section  2,  we  present  some  theoretical  back¬ 
ground.  In  section  3,  we  present  the  algorithm.  In  section  4,  we  prove  results  related  to  the 
convergence  properties  of  the  algorithm.  In  section  S,  we  present  the  updating  scheme  that 
leads  to  a  reduction  in  the  average  number  of  arithmetic  operations  per  step.  In  section  6. 
we  discuss  how  to  initialize  the  algorithm.  Finally,  we  discuss  in  section  7  some 


extensions,  and  in  particular  how  to  extend  our  algorithm  to  solve  convex  quadratic  pro¬ 


gramming  in  O(^fnL)  iterations. 


2.  Theoretical  Background 


In  order  to  facilitate  the  reading  of  this  paper,  we  use  a  notation  roughly  similar  to  the 


one  in  [5].  A  discussion  of  the  main  results  necessary  for  the  development  of  our  algorithm 


is  presented  in  this  section.  These  results  are  adapted  from  [Sj.  A  detailed  discussion  of 


these  results  can  be  found  in  [6]. 


We  consider  the  pair  of  the  standard  form  linear  program  and  its  dual 


( P )  min  cTx 


s.t.  Ax  -  b 


xZO 


(D)  max  bTy 


s.t.  ATy  +  z  =  c 


2  2  0 


where  A  €  R1"*",  c  €  R"  and  b  e  R".  We  impose  the  following  assumptions: 


Assumptions  2.1: 


(a)  The  set  S  ■  |  x  e  R"  ;  Ax  =  b,  x  >  0  [  is  non-empty. 


(b)  The  set  T  ■  j  (y,  r)  e  R"x"  ;  ATy  +  z  ■  c,  z  >  0  [is  non-empty 


(c)  rank  (A)  -  m. 


'.VYVJV'.V.V. 
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Throughout  this  paper,  we  will  denote  a  point  in  R"x  R"x  R"  by  the  lower  case 
letter  w,  i.e., 

w  ■  (x,  y,  z)  e  R"x  Rmx  R" 

The  logarithmic  barrier  function  technique,  usually  employed  in  non-linear  constrained 
optimization,  will  be  applied  to  the  problem  (P).  The  method  consists  of  a  consideration  of 
the  family  of  problems 

ft 

(Pn)  min  cTx  -  n  £  In  Xj 

1 

s.t.  Ax  -  b 
x  >  0 

where  n  >  0  is  the  penalty  parameter.  As  n  -»  0,  we  would  expect  the  optimal  solutions  of 
(P/r)  to  converge  to  an  optimal  solution  of  (P).  This  method,  usually  attributed  to  Frisch 
[1],  recently  came  up  in  [2]  where  an  analogy  with  Karmarkar’s  algorithm  is  presented.  In 
[6],  a  comprehensive  analysis  of  this  approach  is  presented,  where  the  problems  (P)  and  (D) 
interplay  together. 

From  now  on,  a  capital  letter  corresponding  to  a  lower  case  letter  which  denotes  a 

vector,  say  x  *  (X| . x„)T  e  R",  will  denote  the  diagonal  matrix  with  the  components  of 

the  vector  on  the  diagonal,  i.e.,  X  =  diag(x , . x„).  Observe  that  the  objective  function  of 

the  problem  (/V)  is  a  strictly  convex  function.  This  implies  that  the  problem  (Pf. j)  has  at 
most  one  global  minimum,  and  that  this  global  minimum,  if  it  exists,  is  completely  charac¬ 
terized  by  the  Kuhn-Tucker  stationary  condition: 

(i)  ZXt  -  pe  ■  0 

(/i)  Ax  -  b  *  0  ,  x  >  0  (2.1) 

(ill )  ATy  +  z  -  c  *  0 
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where  e  €  R"  denotes  the  vector  of  ones.  In  fact,  we  have  the  following  result: 

Proposition  2.2  :  Under  the  assumptions  (a)  and  (b),  problem  (Pp)  (and  consequently  the 
system  (2.11)  has  a  unique  global  solution  x(p),  for  all  p  >  0. 

Observe  that  in  the  system  (2.1),  x  uniquely  determines  z  from  the  first  equation  and  y 
from  the  third  equation,  since  the  matrix  A  has  full  rank.  For  each  p  >  0,  we  denote  the 
triple  of  solutions  of  (2.1)  by  w(p)  ■  (x(p),  y(p),  z(p)).  Obviously  w(p)  €  S  x  T.  Given 
a  point  w  =  (x,  y,  z)  e  S  x  T ,  we  define  the  duality  gap  at  w  to  be 

g(H-)  =  cTx  -  bTy 

Using  the  two  last  equations  in  (2.1),  we  can  easily  verify  that 

g(w)  =  xrz  ,  w  e  S  x  T  (2.2) 

In  particular,  using  the  first  equation  in  (2.1),  we  obtain  g(w(p))  =  np,  for  all  p  and  there¬ 
fore  g(w(p))  converges  to  zero  as  p  goes  to  zero.  This  implies  that  cTx(p)  and  bTy(p) 
converges  to  the  optimal  value  of  the  problems  (P)  and  (D)  respectively.  In  fact,  we  have 
the  following  stronger  result: 

Proposition  23  :  Under  assumption  2.1,  as  p  -*  0,  x(p)  (  ( y(p ),  z(/i))  )  converges  to  an 
optimal  solution  of  problem  (P)  ((D)). 

The  following  notation  will  be  useful  later.  Let  w  €  S  x  T.  We  denote  by 
/(tv)  *  (/ if**1) . f*(w))T  €  R"  the  n-vector  defined  by 


f,  «  x,z,  ,  i  «  1 . n 

With  this  notation,  the  first  equation  of  (2.1)  written  coordinate-wise  becomes  : 

f,(w(p))  m  x,(p)  z,(p )  =  p  ,  i  =  1 . n  (2.3) 


We  denote  by  T  the  path  of  solutions  w(p)  ,  p  >  0,  i.e.. 


h-00  ■  ?00.  z(p)) ;  p  >0 


The  algorithm  which  will  be  presented  in  the  next  section,  will  "follow"  this  path  r 
with  the  objective  of  approaching  the  desired  solutions  of  the  original  problems  (P)  and  (D). 
The  path  following  procedure  is  described  in  the  next  sections. 

3.  The  Algorithm 

As  in  the  previous  section,  we  denote  a  triple  (jr ,  y,  z)  e  R"x  Rmx  Rn  by  the  lower 
case  letter  w,  i.e.,  w  *  (x,  y,  z).  The  algorithm  will  generate  a  sequence  of  points 
w*  e  5  x  7  ,  k  =  0,1,2,....  where  the  initial  point  w°  is  provided  as  input  to  the  algo¬ 
rithm.  In  this  section,  we  require  that  w°  e  5x7  be  a  point  satisfying  some  criterion  of 
closeness  with  respect  to  the  path  of  solutions  I\  Given  an  LP  problem  in  standard  form,  in 
section  6  '  e  show  how  to  construct  an  equivalent  LP  problem  so  that  assumption  2.1  is 
satisfied.  As  a  consequence  of  this  construction,  we  also  show  bow  to  obtain  an  initial  point 
w°  e  5x7  satisfying  the  criterion  of  closeness. 

Given  a  current  iterate  (x,  y,  z)  €  5  x  7  ,  a  triple  of  directions 
(Ax,  Ay,  Az)  e  R"x  Rmx  R"  needs  to  be  generated  for  the  determination  of  the  next 
iterate.  Throughout  this  paper,  a  triple  of  directions  (Ax,  Ay,  Az)  e  R"x  Rmx  R"  is 
denoted  by  the  symbol  Aw.  Let  (x,  y,  z)  denote  the  next  iterate.  We  obtain  (x,  y,  z)  by 

x  *  x  -  Ax 

y  ■  y  -  Ay 

z  *  z  -  Az 

or  in  more  compact  notation 
w  *  w  -  Aw 

According  to  [5],  the  direction  Aw  chosen  to  generate  the  next  iterate  w,  is  defined  as 
the  Newton’s  direction  associated  with  the  Kubn-Tucker  system  of  equations  (2.1). 
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However,  with  the  objective  of  improving  the  worst-case  complexity  on  the  number  of 
arithmetic  operations,  we  consider  a  slight  variation  of  the  direction  used  in  [5].  If  we 
denote  the  left  uand  side  of  the  system  of  equations  (2.1)  by  H{w)  *  H{x,  y,  z),  the 
Newton’s  direction  Aw  at  w  e  S  x  T  is  defined  by  the  system  of  linear  equations 

DwH(w)  Aw  =  H(w) 

where  Aw  =  (Ax,  Ay,  A z)  e  R"x  R"x  R"  and  DwH(w)  denotes  the  Jacobian  of  H  at 
w  *  (x,  y,  z  ).  We  observe  that  ^(x,  y,  z)  does  not  depend  on  the  argument  y  e  Rm. 
Indeed,  the  Jacobian  of  H  at  w  *  (x,  y,  z)  is  given  by 

Z  0  X 

J(  x,  z)  m  DwH(w )  =  A  0  0 

0  Ar  / 

The  direction  Aw  that  we  are  going  to  consider  is  defined  by  the  following  system  of 
linear  equations 

J(x,  z)  Aw  =  H(x,  y,  z) 

where  the  points  x  e  R"  and  z  e  R"  will  be  chosen  to  approximate  x  e  R"  and 
z  e  R"  respectively  in  a  manner  which  will  be  specified  latter.  More  specifically. 
Aw  =  (Ax,  Ay,  Az)  is  defined  by  the  following  system  of  linear  equations 


Z Ax  +  XAz  =  XZe  -  fie 

(3.1. a) 

A  Ax  =  0 

(3.1 .b) 

A  tA y  +  Az  =  0 

(3.1.0 

where  £  >  0  is  some  prespecified  penalty  parameter.  Observe  that  the  right  hand  sides  of 
(3.1.b)  and  (3.1.C)  are  zero  since  w  is  assumed  to  be  in  the  set  5  x  T.  However  the  right 
hand  side  of  (3.1. a)  is  not  necessarily  zero  and  equals  zero  only  when  the  point  w  lies  on 
the  path  T.  Throughout  this  paper,  the  lower  case  letter  s  will  denote  a  pair 
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(x,  z)  €  R"x  R".  Note  that  the  solution  Aw  =  (Ax,  A y,  Az)  of  the  system  of  equations 
(3.1)  clearly  depends  on  the  current  iterate  w  =  ( x ,  y,  z),  on  the  Jacobian  of  H  at  the 
"approximation"  s  =  (x,  z)  of  s  -  (x,  z)  ,  and  on  the  penalty  parameter  fi  >  0.  In  order  to 
indicate  this  dependence,  we  denote  the  solution  (Ax,  Ay,  Az)  of  the  system  (3.1)  by 

Aw(w,  s,  fi) 

By  simple  calculation,  we  obtain  the  following  expressions  for  Ax,  Ay,  Az 
Ax  =  [z-1  -  Z ~lXA  T(AZ~lXA  t)~'aZ  ~']  (XZe  -  fie) 

Ay  =  -  [(AZ_,jfAr)"lAZ_l]  (XZe  -  fie) 

Az  =  \at  (AZ-xXAt)~'aZ~x ]  (XZe  -  pe) 

Therefore,  to  calculate  the  direction  Aw  s  (Ax,  Ay,  Az),  the  inverse  of  the  matrix 
(AZ~'XAt)  needs  to  be  calculated.  This  is  the  main  motivation  to  consider  just  an  approxi¬ 
mation  s  =  (x,  z )  of  s  =  (x,  z)  so  that  we  do  not  need  to  invert  this  matrix  from  scratch  at 
every  iteration.  If  the  current  diagonal  matrix  Z~xX  differs  from  the  previous  one  by 
exactly  /  diagonal  elements  then,  by  performing  /  rank -one  updates,  we  are  able  to  compute 
the  inverse  of  the  matrix  (AZ~]XAT)  in  0(n2l)  arithmetic  operations.  Observe  that  all  the 
other  operations  involved  in  the  computation  of  Aw  ■  A w(w,  s,  fi)  is  of  the  order  of  0(n2) 
arithmetic  operations. 

We  are  now  ready  to  describe  the  algorithm.  At  the  beginning  of  the  algorithm,  we 
assume  that  an  initial  point  w°  m  (x°,  y°,  z°)  e  S  x  T  is  available  such  that  the  following 
criterion  of  closeness  with  respect  to  the  path  T  is  satisfied: 


II  f(w°)  -  n°e  II  <;  en°  (3.2) 
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where  ||  .  ||  denotes  the  Euclidean  norm,  pn  is  a  positive  constant  and  6  =0.1. 

We  now  state  the  algorithm. 

Algorithm  3.1  : 


Step  0)  Let  w°  €  S  x  T  and  p°  >  0  satisfy  (3.2).  Let  e  be  a  tolerance  for  the  duality 
gap  Let 


6  :=  0.1 


6  :=  0.1 


r  :=  0.1 


Set  k:=0. 


Step  1)  If  cTxk  -  bTyk  <  e,  stop 


Step  2)  Choose  s  =  (J,  :)  in  Rj  x  RJ.  satisfying: 


I  xk  -  x,  I 
I  x,  I 


£  y  ,  i  =  1 . n 


I  a*  ~  I 
I  I 


Sr  ,  i  =  1 . n 


Step  3)  Set  pk*'  :=  /r*(l  -  6  l4n). 

Calculate  Ah’*  ■  Ah’(w'*,  s ,  /t*  +  1). 


Step  4)  Serw*+I  :=  w*  -  Aw*. 

Set  k  :=  k  +  J  and  go  to  step  1. 

In  the  following  sections,  we  prove  that  the  algorithm  above  is  a  valid  one  in  the 
sense  that  it  generates  at  every  iteration  a  point  w*  in  the  set  S  x  T.  We  also  show  that  it 
terminates  in  at  most  0(“^n  max(  loge'1,  log/t,  logjj")  )  iterations.  Finally,  we  present  a 


RRSr 


suitable  choice  for  the  approximation  point  7  =  (7.  7)  (see  step  2  of  the  algorithm)  that 
will  enable  us  to  show  that  algorithm  3.1  solves  the  pair  of  problems  (P)  and  (D)  in  no 
more  than  0(n}  max(  logc'1,  log/i,  log p“)  )  arithmetic  operations. 


4.  Convergence  Results 

We  begin  this  section  by  stating  the  main  result  and  its  consequences.  Given  two  vec¬ 
tors  x  e  R"  and  x  e  R",  we  denote  the  Euclidean  norm  of  the  vector  X-l(  x  -  x  )  by 

II  *  ~  x\\M  •  »  e  . 


‘ 

(  -  ) 

2 

II  X  -  rll,  = 

n  1 

l! 

x  —  x 
•*1 

{  X‘  ; 

The  mam  result  is: 

Theorem  4.1:  Let  w  =  (X,  y.  z)  e  S  x  T  and  p  >  0  satisfy 

|| /(h)  -  pe  ||  5  6p  (4.2) 

Let  7  =  (7,  7 )  €  R"  x  R"  satisfy 
I  x,  -  7,  I 

— - - —  <.  y  .  I  =  1 . n  (4.3) 

I  7,  I 

I  z,  -  7,  I 

— - - -  Sy  •  1  =  1 . n  (4.4) 

I  7,  I 

Let  p  >  0  be  defined  as 

ft  =  p(  1  -  5  /sn  )  (4.5) 

Consider  the  point  w  m  (j,  y,  z)  e  R"x  R"x  R"  defined  by 
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where  Aw  *  Aw(w,  s ,  p).  Then  the  following  holds: 


(a)  The  point  w  is  in  S  x  T  and  satisfies 


II*  -  *11,  *  0.28 


IS  -  *llj  £  0.28 


(b)  |1/(h>)  -  pe\\  £  Op 


(c)  g(w)  *  cTx  -  bTy  £  l.lnp 


The  proof  of  theorem  4.1  will  be  given  at  the  end  of  this  section.  Suppose 


w{p)  a  (x(p),  yip),  z{p))  is  the  point  in  the  path  T  corresponding  to  the  penalty  parameter 


p.  The  criterion  of  closeness  (4.2)  is  then  equivalent  to 


II  /(*’)  ~  /(h’Oi))  ||  £  Op 


since  fiw(p))  =  pe  by  relation  (2.3).  By  theorem  4.1,  relation  (4.9)  will  hold  for  the  new 


point  w  defined  by  (4.6)  and  the  penalty  parameter  p  >  0  given  by  (4.5). 


As  a  consequence  of  theorem  4.1,  we  have  the  following  result: 


Corollary  4.2:  All  points  wk  generated  by  algorithm  3.1  satisfy 


(a)  wk  is  in  S  x  T,  for  all  k  =  1,2,  ...  and 


Jt*+1  -  xk  ||x.  0.28 


z*+l  -  z*  ||,.  <;  0.28 


(b)  ||  f(wk)  -  pke  ||  £  6pk  ,  for  all  it  =  1,2,  ... 


(c)  g(wk)  »  cTxk  -  bTyk  £  1.1  npk  ,  for  all  k  =  1,2, 


where 


=  p°{  1  -  6  hfn  )*  for  k  =  1,2,  ... 
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Proof:  This  result  follows  trivially  by  arguing  inductively  and  using  theorem  4.1.  □ 

By  (b)  of  corollary  4.2,  all  iterates  wk  generated  by  algorithm  3.1  will  satisfy  the  cri¬ 
terion  of  closeness 

II  /(*•*)  -  ||  <;  epk 

Therefore  we  can  view  algorithm  3.1  as  a  path-following  procedure,  i.e.,  the  iterates  h* 
tries  to  trace  the  path  T  so  that  it  eventually  converges  to  an  optimal  solution 
w*  =  (or*.  y*.  z*)  for  the  pair  of  LP  problems  (P)  and  (D). 

We  now  derive  an  upper  bound  on  the  total  number  of  iterations  performed  by  algo¬ 
rithm  3.1. 

Proposition  4_3:  The  total  number  of  iterations  performed  by  algorithm  3.1  is  no  greater 
than  k*  m  [  log(l.lnc_l  p°)  •Jn  I  8  ]  where  e  >  0  denotes  the  tolerance  for  the  duality 
gap  and  p°  is  the  initial  penalty  parameter. 

Proof:  From  (c)  of  corollary  4.2,  we  can  terminate  the  algorithm  as  soon  as 
1.1  npk  £  e 


It  is  enough  to  show  that  k *  satisfies  the  inequality  above.  By  the  definition  of  k% ,  we  have 

log  e  S  -  +  log(l.ln/i") 

sn 


t  k  log 


1  - 


4n 


+  log  (l.l/yO 


=  log 


1.1  np' 


1  - 


4~n 


*  n 


=  log  1.1  np 


k ’ 


Second  inequality  is  due  to  the  fact  that  log  (1  -  x)  £  -x  for  all  x  >  -1  and  the  last  equal¬ 
ity  follows  from  the  definition  of  p  * .  Therefore  k  *  satisfies 


and  this  completes  the  proof  of  the  proposition.  □ 

Let  L  denote  the  number  of  bits  to  represent  the  data  of  the  primal  LP  problem  (P), 
i.e.,  L  is  the  size  of  the  problem  (P).  It  is  a  well-known  fact  that  if,  at  some  iteration  k,  the 
duality  gap  satisfies  cTxk  -  bTyk  £  2-t,  then  we  may  obtain  optimal  solutions  for  prob¬ 
lems  (P)  and  (D)  in  0(mn2)  arithmetic  operations. 

Using  this  observation,  we  obtain 

Corollary  4.4:  If  the  initial  penalty  parameter  p°  satisfies  log  p°  -  0{L)  then  algorithm 
3.1  solves  the  LP  problem  (P)  in  at  most  0(4nL)  iterations. 

Proof:  Follows  directly  from  the  previous  proposition.  □ 

In  section  6,  we  will  see  that  the  initial  penalty  parameter  p°  can  be  chosen  to  satisfy 
log  p°  =  O(L).  One  possible  choice  for  the  approximation  s  m  (x,  z)  on  step  2  of  the 
algorithm  is  to  use  exact  data,  that  is,  to  set  s ,  on  the  *  '*  iteration,  equal  to  sk.  With  this 
choice  of  the  "approximation"  s  ,  we  have  the  following  result: 

Corollary  4.5:  Algorithm  3.1  solves  the  pair  of  LP  problems  (P)  and  (D)  in  no  more  than 
0(n35L)  iterations. 

Proof:  At  every  iteration,  we  need  to  calculate  the  inverse  of  the  matrix  [A(Zk)~'Xk  AT\ 
and  this  requires  0(n3)  arithmetic  operations.  By  corollary  4.4,  algorithm  3.1  terminates  in 

at  most  O(VnL)  iterations.  These  two  observations  immediately  concludes  the  proof  of  the 
corollary.  □ 

In  the  next  section,  we  present  an  updating  scheme  for  the  approximation  s  that  will 
reduce  the  worst-case  complexity  of  algorithm  3.1  to  0(n3L)  arithmetic  operations. 

We  now  concentrate  our  effort  towards  proving  theorem  4.1.  We  should  point  out  that 
some  of  the  arguments  below  become  simpler  when  we  use  exact  data  for  the 
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approximation  s. 

Let  w  =  (x,  y,  z)  e  5  x  T,  7  =  ( x ,  z  )  e  R"x  R"  and  p  >  0.  Let 
Ah’  =  (Ax,  Ay,  Az)  be  the  direction  Ah’(h\  s,  p).  Consider  the  point  w’  defined  as  in  rela¬ 
tion  (4.6).  The  next  result  provides  expressions  for  the  product  of  complementary  variables 
fi(w) ,  i  =  1 . n  and  the  duality  gap  g(w)  =  crx  -  bTy. 

Proposition  4.6:  Let  w  ,7  and  w  be  as  above.  Then  the  following  expressions  hold: 

fi(w)  =  p  +  Ax,Az,  +  (x,  -  x,)Az,  +  (z,  -  z,)Ax,  (4.10) 

(Ax)r(Az)  =  0  (4.11) 

g( H’)  =  X  /.•(*)  =  np  +  (z  -  z)rAx  +  (x  -  x)rAz  (4.12) 

i=i 

Proof:  By  definition,  we  have  for  all  i  =  1 . n 

/,(w)  ■  x,  z, 

=  (x,  -  Ax,)(z,  -  Az(  ) 

=  x(z,  -  (x,Az,  +  z.Ax.)  +  Ax,Az, 

=  x.z,  -  (x,Az,  +  z.Ax,)  +  Ax,Az,  +  (x,  -  x,)Az(  +  (z,  -  z,)Ax,  (4.13) 

Writing  (3.1.a)  coordinate-wise,  we  obtain  for  all  i  =  l,...,n 

(x,Az(  +  z.Ax,)  =  x,z,  -  p  (4.14) 

Relations  (4.13)  and  (4.14)  immediately  imply  (4.10).  Now  we  prove  (4.11).  Multiplying 
expression  (3.1.c)  on  the  left  by  (Ax)r,  we  obtain 

(AAx)rAy  +  (Ax)rAz  =  0  (4.15) 

Relations  (3.1.b)  and  (4.15)  immediately  imply  (4.11).  Note  that  the  first  equality  in  (4.12) 
follows  immediately  from  expression  (2.2)  and  the  definition  of  /.(h).  Summing  expres- 


sion  (4.10)  over  all  indices  i  *  l,...,n  and  noting  (4.11),  we  immediately  obtain  the  second 
equality  in  (4.12).  This  completes  the  proof  of  the  proposition.  □ 

Observe  that  if  x  -  x  and  2=2,  that  is,  the  approximation  is  exact,  then  the  duality 
gap  at  w  becomes  g(w)  =  nfi.  In  this  case,  (c)  of  theorem  4.1  follows  trivially. 

We  now  state  and  prove  some  preliminary  results  that  will  be  useful  in  the  proof  of 
Theorem  4.1. 

Let  w  =  (x,  y,  z)  e  S  x  T,  s  =  (x,  2)  e  R"x  R"  and  (L  >  0.  Let 
Aw  =  (Ax,  Ay,  A 2)  be  the  direction  An^w,  s ,  fi).  We  denote  by  A f  m  (A/,.  ...  ,A f„)T  the 
n-vector  defined  as 


A /  =  (A*,A2, . A x„Azn)T 


(4.16) 


where  Ajt,  and  Ar,  denotes  the  i'h  coordinate  of  the  vectors  Ax  and  A z  respectively.  The 
next  results  provides  an  upper  bound  on  the  Euclidean  norm  of  the  vector  A/. 

Lemma  4.7:  Let  Af  be  defined  as  in  (4.16).  Then,  we  have 


II  A/  ||  S 


(4.17) 


where 


/ mu.  *  min  I  ^2,  I  ;  i  =  1 . n 


(4.18) 


Furthermore,  we  have 


||DA2||2  <,  wMJ—P 

f  min 


(4.19) 


l|5-W  £ 


(4.20) 


where  D  is  the  diagonal  matrix  defined  by 


D  =  (Z-'X)* 


(4.21) 


Proof:  By  equation  (3.1. a),  we  have 


D~'Ax  +  DAz  =  (XZ)*(XZ  -  fie) 


(4.22) 


From  (4.11)  it  follows  that 


(D~'Ax)t(DAz)  =  0 


(4.23) 


Using  the  Pythagorean  theorem,  relations  (4.22),  (4.23)  and  the  definition  of  the  Euclidean 
norm,  we  obtain 

||  5  ~'Ax  ||2  +  ||  DA:  ||2  =  ||  (XZ  )\XZe  -  fie)  ||2 


_  y  C/,-(w>)  -  fi)2 

1*1  i  X  jZ ,  I 


||/(w)  -  fie  II2 


(4.24) 


Inequalities  (4.19)  and  (4.20)  follow  immediately  from  (4.24).  Also  (4.24)  implies  that 


D~'Ax  ||  ||  DAz  ||  £  llf-fol— i*lJl 

2  fm,n 


(4.25) 


On  the  other  hand,  using  the  Cauchy-Schwarz  inequality,  we  obtain 


£  I  Ax,Az,  I  =  £  I  D-'Ax,  I  I  D,jAzj  I 

i*i  i=i 


S  ||  D~'Ax  ||  ||  DAz 


(4.26) 


Since  ||  Af  ||  S  £  I  AjTjAz,  I,  relations  (4.25)  and  (4.26)  imply  the  inequality  (4.17).  This 
i*i 


completes  the  proof  of  the  lemma.  □ 


The  next  lemma  provides  important  relations  that  will  be  useful  in  the  proof  of 
theorem  4.1.  We  observe  that  the  specialization  of  the  results  below  to  the  case  when  the 
approximation  becomes  exact  simply  involves  substitution  of  x  by  x,  :  by  :  and  y  by  0. 

Lemma  4.8:  Let  w  =  (jt,  y,  2)  e  S  x  T  ,  s  =  (x,  z )  e  R"  x  R"  ,  p  >  0  and  ft  >  0  be 
as  in  the  statement  of  theorem  4.1.  Let  A w  ■  (Ax,  Ay,  A:)  be  the  direction  Aw(w,  s ,  ft). 
Let  p  and  q  be  defined  by 


p  =  (1  -  »)(1  +  y)' 


(4.27) 


q  =  (l  +  0)0  -  y)* 


(4.28) 


Then  the  following  relations  hold: 


PP  £  x,z,  <.  qp  ,  1  =  1 . n 


(4.29) 


II  A/  ||  S 


(0  +  S)2l 


(4.30) 


D-’Ax  II2  S  (fl 


(4.31) 


DAz  ||2  £ 


(4.32) 


where  D  is  defined  by  (4.21). 

Proof:  From  (4.3)  and  (4.4)  it  follows  that 


0  <  1  -  y  £  —  £  1  +  y  ,  i  =  1 . n 


(4.33) 


0  <  1  -  y  £  —  £  1  +  y  ,  1  =  1 . n 


(4.34) 


which  implies. 


(1  -  y)2  S  ~r  S,  (1  +  y)2  ,  i  «  1 . n 


v> 
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or  equivalently,  for  all  i«l,...,n 

(1  +  y)~2x,:,  *  x,z,  £  (1  -  y)'2x,z,  (4.35) 

Using  (4.2),  we  obtain 

(1  -  $)n  Z  x*:,  Z  (1  +  0)n  ,i  =  l.  ...  ,n  (4.36) 

Relations  (4.35)  and  (4.36)  imply  that  for  all  i*l,...,n 
(1  +  yr2(l  -  6)n  <.  x,z,  £  (1  —  y)~\  1  + 
which  is  exactly  (4.29). 

Since  ||e||  =  Vn,  relations  (4.2)  and  (4.5)  imply  that 

Il/(w)  -  fze  II2  £  (||  /(w)  -  tie  II  +  II  ne  -  j ie  ||j 

*  (ty  +  \n-fii  lle||) " 

S  (0/i  +  5m)2 

£  (0  +  S)2n2  (4.37) 

Using  lemma  4.7,  relations  (4.29)  and  (4.37)  ,  we  immediately  obtain  (4.30),  (4.31)  and 
(4.32).  This  completes  the  proof  of  the  lemma.  □ 

We  are  now  ready  to  prove  theorem  4.1. 

Proof  of  theorem  4.1: 

(a)  From  (3.1.b),  (3.1.c)  and  the  fact  that  we  S  x  T,  it  follows  that  w  ■  (Jr.  y,  r)  satisfies 
A  x  -  b 
AT y  +  2  ■  c 

We  have  just  to  show  that  x  and  :  are  strictly  positive  to  conclude  that  we  S  x  T.  We 
note  that  (4.7)  and  (4.8)  immediately  imply  that  x  >  0  and  5  >  0  since  x  >  0  and  :  >  0  by 
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7T~T — «* .  I  . 


assumption.  We  will  now  show  that  (4.7)  and  (4.8)  hold.  By  the  definition  of  Jr,  we  obtain 


it  At 

Hi  -  *\\l  =  I  — 


i *  I  l  **i 


(4.38) 


where  D  =  diag(Du . Dnn)  is  given  by  (4.21),  i.e., 


D„  =  —  i  =  1 . n 


(4.39) 


Using  (4.2),  (4.3),  (4.4)  and  the  expression  for  £>()  above,  we  obtain  for  all  i  =  1 . n 


A,  *i 


UJ  U*.- 


1  -  y  (1  -  e)n 


The  last  inequality  together  with  (4.38)  implies 


(4.40) 


iiJr  -  jtii2  s  *  y  - - - 

x  l  -  y  (l  -  0)M 


D~]Ax  H2 


(4.41) 


Inequality  (4.31)  of  lemma  4.8  and  (4.41)  then  imply 


»  -  s  <, <e  +  i) 


(4.42) 


Substituting  in  (4.42)  the  values  of  6,  S,  y,  and  p  from  (3.3),  (4.27)  and  (4.28)  we 
(4.7).  In  a  similar  way,  one  can  prove  (4.8).  This  completes  the  proof  of  (a). 

(b)  Writing  expression  (4.10)  in  vector  notation  and  noting  (4.16),  we  obtain 
f(w)  -  fie  =  Af  +  (X  -  X)A :  +  (Z  -  Z) Ax 
By  the  properties  of  norms,  we  obtain 

||  f(w)  -fie  ||  $  II  Af  ||  +  ||  (X  -  X)A:  ||  +  ||  (Z  -  Z)Ax  || 


Using  the  definition  of  the  Euclidean  norm,  we  obtain 


(X  -  X)Az  ||2  =  £  (x,  -  x,)\A2,)2 


i=i 


=  £  [Du\x,  -  Xi)UDu A:,)2 


i= i 


where  Du  is  as  in  (4.39).  Using  (4.3),  (4.29)  and  (4.39),  we  obtain 


[Dti '(?,  -  *,)]  =  (?,?,) 


(  \ 5 
-  x, 


Z  qpy2  ,  i  =  l . n 


In  a  similar  way,  one  can  prove  that 


I Du(?i  -  z,)]2  Z  qpy2  ,  t  =  1 . n 


Relation  (4.32)  of  lemma  4.8  and  relations  (4.44)  and  (4.45)  imply  that 


(X  -  X)Az  r  * 


q-f{e  +  5)V 


obtain 


(4.43) 


(4.44) 


(4.45) 


(4.46) 


(4.47) 


In  a  similar  way,  one  can  prove  that 
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From  (4.43),  (4.30),  (4.47)  and  (4.48)  it  follows  that 

||  f(w)  -  fie  ||  £  +  2*0  +  5)(-2)*  p 

2 P  P 

From  the  expression  of  p  in  (4.5),  we  obtain 

^  (1-5) 

which  implies 

||  f(w)  -pe\\<,  — ~  ^  S)--  +  2*0  +  5)(  V  p 

(1  -  6)  [  2p  p 

Substituting  the  values  of  0,  8,  y,  p,  and  q  from  (3.3),  (4.27)  and  (4.28)  we  obtain 

||  /(*')  -  pe  ||  £  Op  (4.49) 

and  this  completes  the  proof  of  (b). 

(c)  The  first  equality  of  expression  (4.12)  says  that 

£(*’)  =  eTf(w) 

Using  the  Cauchy-Schwarz  inequality  and  expression  (4.49),  we  obtain 

g(  *>  s  Ikll  II/(h)II 

£  Jn  (  II /(H’)  ~  M^ll  +  \\pe\\  ) 

J»  *fn  (Op  +  Jnp) 

£  (1  +  0)np 
-  \.lnp 

since  0  =  0.1.  This  completes  the  proof  of  (c).  □ 
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5.  A  Good  Choice  for  x  and  z 

In  this  section,  we  show  how  to  choose  the  approximation  s  *  (x,z)  e  R"x  R"  of 
sk  =  (jc* ,  zk)  in  step  2  of  algorithm  3.1  so  that  its  worst-case  complexity  reduces  to  the 
order  of  0(n3L)  arithmetic  operations. 

The  choice  of  the  approximation  s  is  made  by  an  updating  scheme  as  follows  (In  the 
procedure  below,  k  stands  for  the  iteration  count): 

Updating  scheme  5.1: 

For  k  :=  0,  set  x  :=  x°  and  z  :=  z° 

For  k  >  0  do 

For  i  =  1 . n  do 

If  one  of  the  following  holds 


(a) 


(b) 


I  xk  -  x,  I 

I  X,  I 


>  y 


i  ' 

I  Z.  I 


I  >  y 


then  set  x,  :=  xk  and  z,  :=  zk. 

In  section  3,  we  saw  that  the  main  effort  in  each  iteration  was  in  computing  the 
inverse  of  the  the  matrix  (AZ_lXAr).  If  only  /  diagonal  elements  of  the  matrix  {Z''x ) 
change,  this  computation  can  be  carried  out  in  0(n2l)  arithmetic  operations  by  means  of  / 
rank -one  updates.  If  we  use  the  updating  scheme  5.1  then  the  ith  diagonal  element  of  the 

matrix  (Z~'X)  changes  only  when  inequality  (a)  or  (b)  of  scheme  5.1  is  satisfied.  Next,  we 
provide  an  upper  bound  on  the  total  number  of  diagonal  element  changes  that  occurs  on  the 

matrix  (Z~xX)  during  K  iterations  of  algorithm  3.1.  Consider  the  following  two  sets: 


S  r  =  *  ;  - >  y .  1  *  *  <.  K 

I  x,  I 


Tr  -  i  k 


I  2*  -  2, 


I  2.  I 


>  y  ,  \  <  k  <  K 


where  K  is  the  number  of  iterations  performed.  Note  that  by  corollary  4.2.  the  iterates 


h*  =  (xk,  y*.  2*)  generated  by  algorithm  3.1  satisfy 


v*+l  -  x*||,.  S  0.28 


|r*+l  -  2*11..  <  0.28 


for  all  integers  k  2  0. 


The  following  result  can  be  proved  by  slightly  modifying  the  arguments  in  section  5 


of  (3). 


Proposition  5.2:  Let  S  and  T  be  defined  as 


s  =  £is*i 
1  =  1 


T  =  £l  7,*  I 


Then,  we  have 


S  £  4.5 -fnK 


T  £  4.5  JnK 


As  a  consequence  of  this  result,  we  have  : 


Corollary  5.3:  Algorithm  3.1  coupled  with  the  updating  scheme  5.1  solves  problems  (Pi 


and  (D)  in  no  more  than  0(n*L)  arithmetic  operations. 


Proof:  From  corollary  4.4,  we  know  that  algorithm  3.1  finds  an  optimal  solution  in  O(-JnL) 
iterations.  Proposition  5.2  implies  that  the  total  number  of  rank -one  updates  is  then  of  the 
order  of  O(nL).  Since  each  rank-one  update  involves  0(n arithmetic  operations,  the  total 
number  of  arithmetic  operations  is  then  of  the  order  of  0(nyL).  This  completes  the  proof  of 
the  corollary.  C 

6.  Initialization  of  the  Algorithm 

Given  an  LP  in  standard  form  and  its  dual,  we  have  assumed  in  section  2  that  condi¬ 
tions  (a)  and  (b>  of  assumption  2.1  were  valid.  In  general,  this  is  not  the  case.  In  this  sec¬ 
tion,  we  show  how  to  transform  an  LP  problem  in  standard  form 

(P )  min  rrv 

s  t  Ax  -  b  ,  A  full  row  rank 

v  2  0 

to  an  equivalent  LP  problem  in  standard  form  so  that  the  new  problem  and  its  dual  satisfy 
the  assumption  2.1.  As  a  result  of  the  process  to  be  described  below  ,  we  will  also  obtain  an 
initial  point  satisfying  the  criterion  of  closeness  (3.2).  The  algorithm  can  then  start  from  this 
initial  point. 

Let  m  and  n  be  integers  such  that  A  e  First  we  scale  the  problem  ( P ) 

by  introducing  a  change  of  variables  v  =  AJ  where  A  >  0  is  a  scaling  factor  and 
x  e  R'""2’.  The  scaling  factor  A  is  chosen  large  enough  so  that  any  basic  feasible  solu¬ 
tion  v  of  (P)  satisfies  e T v  <,  (n- 2)A.  We  then  obtain  the  following  equivalent  problem 
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min 


(XcT)  x 


s.t.  (XA  )  x  =  b 

eTx  +  7n_,  =  n-1 
x  £  0  ,  ;>  0 

where  xn_i  is  a  slack  variable. 

Now  we  want  the  vector  of  ones  to  be  a  feasible  solution  of  the  new  problem.  We 
introduce  an  artificial  variable  x„  with  a  large  cost  M  to  obtain 

min  ( XcT )  x  +  M  x„ 

s.t.  (XA )  x  +  (b-XAe)  x„  -  b 

eTx  +  +  x„  =  n 

x  £  0 , 3c„_,  £  0  ,  x„  £  0 

Letting  B  =  [XA  I  b  -  XAe  I  OJ  e  R('"*l)x",  cT  =  (XcT ,  0,  M)  and 
x  -  (J,  xn)  e  R",  we  can  rewrite  the  previous  problem  as 

( P )  min  cTx 

s.t.  Bx  =  b 


T 

ex  =  n 


x  £  0 


a 


We  can  recast  this  problem  in  the  notation  of  problem  (P)  of  section  2  just  by  letting 
A  c  Rmxn  and  b  e  Rm  be 


9 


A  = 


. 

B 

h  — 

b 

eT 

0  — 

n 

\ 
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We  now  show  that  this  problem  satisfies  assumption  2.1  of  section  2.  Note  that,  by  the 
construction  above,  the  vector  of  ones  e  e  R"  is  feasible  for  problem  (P).  This  shows  that 
problem  (P)  satisfies  (a)  of  assumption  2.1.  By  assumption,  the  matrix  A  has  full  row  rank. 
This  immediately  implies  that  the  matrix  B  has  hill  row  rank.  Since  Be  -  0,  it  follows  that 
the  matrix  A  has  full  row  rank.  Therefore,  problem  (P)  satisfies  (c)  of  assumption  2.1.  To 
see  that  problem  (P)  also  satisfies  condition  (b)  of  assumption  2.1,  consider  its  dual 

(D)  max  bTri  +  nt 

s.t.  BTr\  +  et  +  2  =  c 
z  Z  0 

where  rj  e  R('"_l),  /  e  R  and  the  slack  z  c  R".  Problem  (D)  has  a  feasible  point  with 
the  corresponding  slack  z  strictly  positive.  To  see  this,  set  r)  =  0  and  choose  any  r  satisfy¬ 
ing 

r  <  min  c, 

< 

Hence  condition  (b)  of  assumption  2.1  is  satisfied.  We  will  now  obtain  an  initial  point  satis¬ 
fying  the  criterion  of  closeness  (3.2).  Let  x°  =  e,  n  0  =  0,  and  t°  =  -  n°  where  n°  satisfies 

(6.1) 

Then  x"z^  =  c;  +  n°  ,  j= 1 . n  and  the  criterion  for  closeness  (3.2)  becomes 

(  n  * 

£  cf  SO  n° 

which  is  satisfied  due  to  (6.1).  In  summary,  we  can  apply  the  algorithm  of  section  3  to 
solve  the  pair  of  problems  (P)  and  (D)  with  the  initial  point  w°  *  ( x° ,  y°,  z°)  determined  by 
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y°  m  (W°,  t°)  =  (0.  -n°) 
z°  =  c  +  n°e 
where  n°  satisfies  (6.1). 

7.  Remarks  and  Extensions 

The  following  observations  are  in  order: 

(1)  The  purpose  of  this  paper  is  to  present  a  theoretical  result.  Thus  in  order  to  simplify  the 

presentation,  we  constructed  fi  =  p(l  -  6  iJn).  Obviously,  one  can  use  ft  which  is  less 
than  or  equal  than  the  above  one,  but  still  satisfying  (b)  of  theorem  4.1  and  relations  (4.7) 
and  (4.8).  hi  this  way,  one  can  accelerate  the  convergence  of  the  algorithm. 

(2)  Additional  improvements  in  actual  implementation,  which  are  possible,  such  as  more 
judicious  selection  of  0,5  and  y,  together  with  actual  test  results,  are  the  subject  of  a  forth¬ 
coming  paper. 

7.1.  Extensions  to  Convex  Quadratic  Programming 

The  same  ideas  presented  here  can  be  applied  to  the  convex  quadratic  programming  as 
follows.  Let 

( Q )  min  cT x  +  ~  xT Qx 

s  t  Ax  =  b 
x  2  0 

where  Q  c  R"x"  is  positive  semidehnite.  Using  the  same  approach  and  following  [6]  one 
can  define  the  following  barrier  problem 


with  the  Kuhn-Tucker  conditions: 

(/  )  ZXe  -  ne  =  0 

(ii)  Ax-b=0,x>0 

(i'll )  A  ry  +  z-  c- £)*  =  () 

As  in  the  linear  case  the  path  of  solutions  to  ( Qn )  converges  to  the  optimal  solution  of  (Q) 
as  n  tends  to  zero  (see  Meggido  [6]).  Similarly  to  algorithm  3.1,  given  a  point 
h’  =  (jr,  y,  z),  one  can  construct  Aw  =  (Ax,  Ay,  A z)  by  solving  the  following  system  of 
linear  equations 

ZAx  +  XAz  =  XZe  -  ft  e 

A  Ax  =  0 

At Ay  +  A:  -  QAx  =  0 

and  then  proceed  to  w  by 
w  =  w  -  Aw 

It  can  be  shown  that  convergence  properties  of  this  algorithm  are  the  same  as  the  ones 

of  algorithm  3.1  (  i.e.,  0(4nL)  iterations,  each  of  which  requires  at  most  0(n 3)  arithmetic 
operations).  The  proofs  are  slight  modifications  of  the  proofs  of  this  paper.  The  details  of 
this  algorithm  together  with  implementation  tests  are  the  subject  of  a  forthcoming  paper. 
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